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ABSTRACT 


ROBERT S. BERNARD, Doctor of Philosophy, 1981 

Major: Engineering, Department of Aerospace Engineering 

Title of Dissertation: Approximate Factorization for 

Incompressible Flow 

Directed by; Dr. Joe F. Thompson 

Pages in Dissertation: 92 Words in Abstract; 250 

Abstract 

The technique of approxlroate factorization (AF) , as forismlated 
by Beam and Warming, is employed for computational solution of the 
incompressible Navier -Stokes equations. In each time step, the AF 
algorithm is used to solve the vectorized momentum equation in delta 
form, based on the pressure calculated in the previous time step. The 
newly-calculated velocities are substituted into the pressure equation 
(obtained from a linear combination of the continuity and momentum 
equations) , which is then solved by means of line SOR. 

The combined AF/SOR algorithm is developed for arbitrary curvi- 
linear coordinates in two dimensions, but can easily be extended to 
three. The coordinates used for computation are the body-fitted coor- 
dinates of Thompson et al. , generated numerically by the more recent 
method of Thompson and Mastin. 

Computational results are presented for the NACA 66^018 airfoil 
at Reynolds numbers of 1000 and 40,000 and attack angles of 0 and 6 de- 
grees. Comparison with wind-tunnel data for Re - 40,000 ijidicates 
good qualitative agreement between measured and calculated pressure 

distributions. Quantitative agreement is only fair, however, with the 
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calculations somewhat displaced from the measurements. Furthermore, 
the computed velocity profiles are unrealistically thick around the 
airfoil, due to the excessive amount of artificial viscosity needed 
for stability. 

Based on the performance of the algoritlm with regard to stability, 
it is concluded that AF/SOR is suitable for calcvlations at Reynolds 
numbers less than 10,000. Speedwise, the method is faster than point 
SOR by at least a factor of two. 
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CHAPTER I 


INTRODUCTION 

Since the late 1960'e, computational fluid dynamics (CFDJ has 
evolved as a distinct branch of fluid mechanics, due to numerous advan- 
ces in computer technology and numerical methodology. Now, with the 
increasing cost of fuel and hardware and the decreasing cost of elec- 
tronic calculation, CFD offers a viable alternative to laboratory 
tes ting . Chapman [1] gives a concise overview of the state-of-the-art 
(as of 1979), summarizing its past development as well as its progno- 
sis for the future. 

Computational solutions of the Reynolds-averaged Navier-Stokes 
equations date back to the early 1970' s. At first the machine- time re- 
quirement was prohibitive, and calculations were restricted to simple 
geometries and low Reynolds numbers. These initial limitations were 
greatly reduced, however, with the development of improved algorithms 
such as MacCormack's rapid solver [2] and effective grid-generation 
techniques such as the body-fitted coordinates of Thompson et al. [3]. 
Further reduction of the machine-time requirement was achieved when 
Beam and Warming [A] and Briley and McDonald [5] independently formu- 
lated the technique of approximate factorization (AF) . 

The AF algorithm is an alternating-direction implicit (ADI) scheme, 
suitable for partial differential equations of the parabolic/hyperbolic 
type. Although it requires more computer storage than other methods, 
it is also much faster, maintaining second-order accuracy with one iter- 
ation per time step. Furthermore, since it is fully implicit, AF is not 
subject to the restrictive time-step limitation that impedes explicit 





methods. 


To date, approximate factorization has been successfully employed 
in a number of applications Involvtog supersonic, transonic, and com-’ 
pressible subsonic flow [6, 7, 8], Except for the work of Steger and 
Kutler [9], however, there seems to have been no published calculation 
for incompressible flow. This is not surprising, since the incompres- 
sible equations of motion are not amenable to the AF algorithm, at 
least not without some modification. 

The spatial factorization used in the ADI scheme requires that 
each of the governing equations contain a time derivative. Unfortu- 
nately, the time derivative in the continuity equation is lost in the 
Incompressible limit. Steger and Kutler circumvent this difficulty 
artificially by addiiig a pressure time derivative to the continuity 
equation: 

+ 8’^‘u - 0 (1.1) 

<3t 

whei*e 8>>1 so that 7>u -’>■ 0 as the flow approaches the steady state. 

The calculated transient flow is, however , only asymptotically correct, 
since Eq, (1.1) is a fabrication and not a true equation of motion. 

Available methods for calculating transient incompressible flow 
seem expensive and inefficient compared to those for the compressible 
case, especially considering that the incompressible equations are 
"simpler". One is forced to choose between explicit techniques , lim- 
ited to small time steps, and implicit techniques, requiring many iter- 
ations . The very existence of an efficient algorithm such as approxi- 
mate factorization prompts one to think that an equally efficient recipe 

2 










might be possible for incompressible flow» That, at least, was the 
motivation for the investigation reported herein. 

The AF technique will not accomodate the incompressible continuity 
equation without a time derivative. Nevertheless, it is quite appli- 
cable for the incompressible momentum equation as long as the pressure 
is specified by some other means. The implementation of separate al- 
gorithms for the momentum and continuity equations allows one to take 
advantage of approximate factorization for the velocity calculation but 
still requires some sort of iterative solver for the pressure. This 
also precludes the simultaneous calculation of velocity and pressure, 
so that the combined method is not fully implicit. 

In order to investigate the utility of a two-phase solver based on 
approximate factorization, the two-dimensional incompressible momentum 

4t 

equation^ in conservative form is first transfonned from cartesian co- 
ordinates to arbitrary curvilinear caordinates. The transformed 
equation is discretized by means of a three-point backward difference 
in time and central difference in space, so that the difference oper- 
ator can be factored into a product of two unidirectional operators. 

The factored difference equation is then solved in an ADI sequence with 
the aid of a direct inversion for each of the two coordinate directions. 
Throughout this entire process, the pressure is assumed to have been 
previously calculated. 

A Poisson equation for the pressure is obtained from a linear com- 
bination of the continuity and momentum equations. After discretization 

t The terms ''momentum equation" and "Navier-Stokes equations" are inter- 
changeable for Incompressible flow. 
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with central difference approximations > successive over-relaxation 
by lines (line SOR) is used to solve the pressure equation » Line 
SOR was chosen as the pressure algorithm because it facilitates the 
implicit treatment of Neumann boundary conditions. 

The body-fitted coordinate system is generated numerically for 
a single-eleiuent airfoil using the method described in Thompson 
and Mastin [lO]. Subsequent calculations for the velocity and pres- 
sure are executed in the transformed plane, rather than the physical 
plane* on a 113 x 51 rectangular grid with unit mesh spacing. Due to 
the one-to-one mapping between the physical and transformed planes, the 
recovery and presentation of results in the physical plane is straight- 
forward, 

A linear gradual start, equivalent to a uniform body force or stat- 
ic pressure gradient, is used to accelerate the flow from rest to its 
ultimate freestream velocity. At each time step the (explicit) pres- 
sure solver is executed first, based on velocities from the previous 
time step. Then the (implicit) AF momentum solver is executed, gener- 
ating velocity increments for the current time step. 

An algebraic turbulence model, formulated by Baldwin and Lomax 
[14], is employed for calculating eddy viscosity in the Reynolds-averaged 
momentum equation. In addition, an artificial viscosity, proportional 
to the divergence of the velocity, is appended to the physical viscos- 
ity in order to help dissipate numerical oscillations. 

The foregoing paragraphs summarize the form and content of the 
AF/SOR method. The questions to be addressed are:' 

4 




1. Is this inetUod stable and aecurate for the same Reynolds 
numbers and time increments attainable vith conventional 
implicit methods such as point SOR? 

2. What are its limitations in terms of high Reynolds 
number » time-step size, and boundary conditions? 

3. Under what circumstances is the method preferable to 
conventional Implicit methods? 

The answers are given in Chepter IX, 





CHAPTER It 


EQUATIONS OF MOTION IN CARTESIAN COORDINATES 
2,1. Navler-Stokes Equations 

Newton's Second Inw (conservation of momentum), applied to the 
fluid Inside an arbitrary material volume V, takes the mathematical 
form * 


Dt 


(pu^)dV + / (pu^)u^n^dS 


( 2 . 1 ) 


S J V 

Equation (2.1) represents the well-known Navler-Stokes equations In 
their integral form. The index "1" denotes any of the three cartesian 
directions, X|,X 2 ,X 2 » end the Einsteiii summation convention has been 
employed for the index "j", The dimensional variables are 

p <= density 
- velocity 

S = material surface 

n^ - unit vector, normal to S 

5 . , = shear stress tensor 
ij 

p = pressure 

6.. = Kroneoker delta 

ij 

= body-force acceleration, due to 
gravity, non-inertial coordinate 
system, etc. 

The divergence theorem transforms Eq. (2.1) to 




( 2 . 2 ) 


/ ■^(puj^)dv + / ^(pu^Uj)dV 


/. (1 


_ lE. 


+ P8j)dV 


V id i 

From a computational stond|)Olnt» the integral relation (2.2) is prob- 
ably the most important of the many possible forms of the Navier-Stokes 
equations. In order to achieve a meaningful representation of Newton’s 
second law, numerical simulations should begin with Eq. (2.2), pre*- 
servlng the integral fornv as accurately as possible. 

Since the volume V is arbitrary, the integral operators con be 
removed from Eq. (2.2), leaving only the relation between the inte- 
grands, which is 

3o 


-it + jr ^ hs 

J J 




(2,3) 


For an incompressible fluid the density is constant, so Eq. (2.3) be- 
comes 

a „ 

9P 


3 


p '' 3x^ ^ij 3x^^ ^ % 


and Stokes’s Law gives the shear stress as 

3u, 3u 

^13 “ 




(2.4) 


(2.5) 


where y is the viscosity of the fluid. 

At this point it is expedient to nondimensionalize Eqs. (2,4) and 
(2.5) . Accordingly, the dimensional variables are now replaced by 
their noodimenslonal counterparts: 

(2.b) 




Uj /U 

i “ 
x^/i 

tu /i 

icq* 


(2.7) 

( 2 . 8 ) 
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m/m, 




(2.9) 

( 2 ) 10 ) 

( 2 . 11 ) 


p " p/(puj) 

h " 

The reference quantities are 

" freestream velocity 

f - characteristic length (in this case, airfoil 

chprcl) 

M„ " freestreaitt viscosity 

Gombining Eqs. (2,4) and (2,5) with (2.6) - (2.11) the dimensionless 
Navler-Stokes equations for incompressible flow are 


at ax. ^ J Dx. ^ ax. ax. ax, 

j 1 j j 1 

where Be is the Reynolds number, given by 


( 2 . 12 ) 


Re 


M„ 


(2.13) 


The viscosity is a material constant for Incompressible, noncon- 
ducting fluids, It is retained liere as a variable, however, in order 
to facilitate the implementation of an algebraic model for turbulence. 
2,2 Continuity Equation 

In addition to the Navler-Stokes equations for conservation of mo- 
mentum, all fluids must obey the continuity equation, which expresses 
conservation of mass: 


/ 1^ dV + / pu.n.dS 0 
V s J a 


(2.14) 


Applying the divergence theorem and eliminating the volume integrals as 
before, Eq. (2,14) reduces to 




(2. IS) 




which, for incompressible fluids, is simply 

« 0 (2,16) 

j 

Equation (2,16) is unchanged by the introduction of nondimenslonal 
variables , 

2,3 Nondimensional Equations of Motion in Two Dimensions 

From this point on, all variables used will be nondimensional, and 
the circumflex (*) will be dropped from the notation. Identifying x, 
y, u, V with Xj^, X 2 , u^, U 2 respectively, the equations of motion for 
incompressible flow in two dimensions are 

+ (p + u^) + (uv) = —[2(pu ) + (mu ) + (mv ) ] (2.17) 

c X y A X y y ^ y 

V + (uv) + (p + v^).. « -^[(mu ) + (mv^) + 2(mv ) ] (2.18) 

t X yKeyx Xx 7y 

u + V ■= 0 (2.1.9) 

X y 

Equations (2.17) and (2.18) represent the x- and y- components of the 
momentum equation (2.12), and Eq. (2. 19) replaces the continuity 
equation (2*20). Ihe subscripts x, y, and t indicate x-, y-, and t- 
derivatives. The conservative form (see section 3.2) has been retained 
in (2.17) and (2,18) to preserve the integral relation (2.2) in the 
finite-difference algorithm. 


t In the absence of a body force. 


CHAPTER III 


CURVILINEAR COORDINATES 
3 . 1 Coordinate Generation 

Whenever a physical problem Is to be solved by analytical means » 
one of the first steps Is the choice of an expedient coordinate sys- 
tem - preferably one that will exploit any symmetries and keep the 
boundary conditions simple. If a problem yields an exact solution^ it 
is usually because there exists a coordinate system in which the gov- 
erning equations are algebraically tractable. 

From a computational standpoint, geometry is equally Important. 
Moreover, except for some additional bookkeeping, the formulation of 
difference equations is no more difficult in curvilinear coordinates 
than it is in cartesian coordinates. The practical obstacle is the 
fact that engineering problems seldom involve shapes for which coordi- 
nate systems can be specified analytically. Thus the difficulty is 
mainly one of finding or generating a coordinate system with boundaries 
that match the physical boundaries. 

Thompson et al. [3] have developed a method of grid generation 
whereby coordinate systems can be specified for arbitrary geometries. 
Consider, for example, an airfoil cross-section in two dimensions 
(Figure 1). The desired curvilinear coordinates, C and n , should have 
the following properties ; 

1. The coordinates § and n must be single-valued functions 
of the cartesian coordinates x and y. 

2. Extrema in g and n may occur only on the boundaries, or on 
a branch cut, where the values of both coordinates are 
specified. (In this case n is fixed and C varies along 
the airfoil surface) . 

10 : 






Now it happens that the deslted properties oi 5 and n ath the 
same as those o£ solutions to elliptiR partial differentiaX equations 
with Dirichlet houndary conditions* The simplest of these is Laplace’s 
equation, which allows the values o£ 5 and q to he specified at will 
along the boundaries, A more useful equation, however, is Poisson's 
equation, which facilitates control of both the mesh spacing and the 
cell aspect ratio . Thus lines of constant f and/or n can be densely 
concentrated In regions where large physical gradienta are expected. 

In two diniensions, a system of two Poisson equations must be 
solved to obtain C(x,y) and n(x,y) : 


XX 'y: 

^'‘c 

' C 

(3.1) 

q + q 

XX yj 


(3.2) 


Depending on their magnitude and spatial variation, the attraction 
functions P(5,q) and Q(f,n) can be used to eoncentrate coordinate 
lines about a given point or curve. Since the numerical boundary con- 
ditions can be incorporated more accurately in the rectangular Cn-plane 
(Figure 1), it is expedient to solve the transformed system of 
equations: 




(3.3) 

V« ■ + Vnn ' ' 


(3.4) 


Although Eqs, (3,3) and (3,4) are more complicated than (3.1) and (3,2) , 
computation is far simpler on the rectangular grid. The coefficients 


in Eqs. (3.1) - (3.4) are themselves functions of g and q, given below: 




U 




(3.5) 


■ Vn * 


- *1 + y| 


■'c ■ Vn “ Vs 


(3.6) 

(3.7) 

(3.8) 


Equeiblotts (3 *3) and (3.4) are solved by point SOR after the at« 
traction functions have been specified, fhe procedure for Q io [lO]} 


'<") • 'j + o('S)-H ” - h'’ 


QCCrn) 


h + r* y - 
*- r,^ 


* 6 .r r V. 

j 2 V “ r ' 
c ^ 


(3.9a) 

(3.9b) 

( 3 . 10 ) 


with b chosen so as to produce a specified value of r' (1) . 

The first n"liue (n=l) represents the body, and r^ is one-half chord 


length (the radius of a circle circumscribing the airfoil cross-section) , 
The quantity r. is the radius of a circle tangent to the outer bound- 
ary, n=J. In the solution of Eq. (3.2) , the effect of Q defined by 
Eqs. (3.9) and (3.10) is to position the line q®N at a distance roughly 
equal to r^^ - r^ from the body. For the second n- line (na2), this is 
about 1% of the boundary-layer thickness for the Blasius fiat plate 
solution, i.e. , 

r*(l) = r„ - r, = 0.0l(-“) (3.11) 
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$.2 Coordlnatsa Transformation 

iii j — ■ « ii iL ii » , i» n ■■ u ii r iiiii i iiii ^« «i« i i w w w« « ^i^ ^ 

In Chapter lit the dimensionless Navler-Stokes equations (2fl7) 
and (2»18) are given in physically conservative form. That is, deriv- 
atives such as (u^) are not expanded into the analytically equivolent 
form 

(u2)j^ - 2uu^ (3.12) 

The reason for this decision is thot finite-difference approximations 
for the conservative expressions are consistent with the Integral form 
of the Navier-Stokes equations (2.2), while those for the nonconserva- 
tive expressions are not. Koache demonstrates this at length in his 
hook [ll]. Suffice it here to say that the conservative form is go';i- 
erally safer for computational work. 

A similar consideration arises when the Navier-Stokes equations 
are to be transformed from cartesian coordinates to curvilinear coor- 
dinates. If the chain rule alone is used to expand derivatives such as 


"x ” «x “5 ^ "x “n 


(3.13) 


then the finite-difference equations will be inconsistent with the in- 
tegral Navlcr-Stokes equations. Thus it is Important that the trans- 
formation rule be geometrically conservative in the sense that Eqs. 

(2.17) and (2.18) are physically conservative. The transformation given 
below has the desired property [8]. 

Let Aj^, ^3 the components of the vector A in cartesian co- 
ordinates (xj, Xg, Xg), Employing the Einstein summatiott convention, 
the divergence of A in curvilinear coordinates (|j, takes 

the fosm 








and J denotes the Jacobian 


i J "J 


3(5^,12,53) 

3 (x 3,X2,X3) 


(3. 16) 


In Chapters IV and V the two-dimensional coordinates x^ and are 


Xj = X 


x„ = y 


“ I 




(3.17) 


and the Jacobian is 


= I n - C n 

X y y X 


(3.18) 


I 1 

Applying Eqs. (3.14), (3.15), (3. 17) , and (3. 18) in two dimensions, the 
divergence of A becomes 


9A, 8A„ 3A. 3A„ 

^ X “ ss T ( i. 4 . 

3x 3y ^ 35 9n V 


(3.19) 


where 


^1 “ J^^x ^1 ^y ^2^ 
'' 1 

= T(n.. A, + n. A„) 


(3.20) 


(3.21) 


t The Jacobian J is not to be confused with the maximum n-value (r)=J) 
in section 3.1. 


>— ~ .W, , ii, : -i| f -f .li.j . 



The X- or y- derivative of a scalar is the same as the divergence of 
a vector with only one component; for example, 


/V A 


* J ( Uj, + V ) 

' ? n 

(3.22) 

A £ u 

u - 

J 

(3.23) 

n u 


(3.24) 


Second 'derivatives should be calculated via a combination of the chain 
rule (3.13) and the transformation (3.14), as shown in the example be- 
low 


u « (u ) 

XX X 


(u ) = j[(u ) + (v ) ] 

X X ^ ^ n 


A 5 U 5 

_ ^X X _ ^X/v , V 

"x - "T- - "5 + \ 

- ^ “x "x,, 

''x “ “T“ ■ T<«x "5 \ "n> 


(3.25) 

(3.26) 

(3.27) 

(3.28) 


Note that Eq. (3.13) has been substituted for u^ in Eqs. (3.27) and 
(3.28), since only the outside derivative in (3.26) need employ the 


transformation (3.14) to preserve the Integral relation (2.2). 


CHAPTER IV 


SOLUTION OF THE NAVIER-STOKES EQUATIONS IN CURVILINEAR COORDINATES 
4.1. Transformed Equations 

Using the transformation (3.14) for an arbitrary coordinate sys- 
tem, the X- and y- momentum equations (2.17) and (2,18) can be ex- 
pressed as a single vector equation: 

T I? + w ^ i ^ 

where 

C = ?(x,y) (4.2) 


n = n(x»y) 

j = 5 n - ? n 0 

^x 'y y 'x 


q = 


A(q) = 
B(q) = 


u 

V_ 

’a^(q)"' 

a2(q) 

bj^(q) 



w 
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w 


12 


w 


22 


M2 


X 


22 


'12 


'22 


(4*3) 

(4.4) 

(4.5) 

(4.6) 

(4.7) 

(4.8) 

(4.9) 

(4.10) 


The individual components of the vectors (4 >6) and (4.7) and matrices 
(4.8) - (4.11) are as follows: 


ai(q) “ 

(4.12 

a£(q) = + 5y(p + v^)] 

(4.13 

bj(q) = jCnj^(p + u2) + HyUv] 

(4 . 14 

b2(q) = + hy(p + v2)] 

(4.15 

Wii = i(2?2 + ?2) 

(4.16 

w, o = W-, = yC ? 

12 21 X y 

(4.17 

w„„ = y(?2 + 2 ? 2 ) 
22 y 

(4.18 

“u ■ ''ll “ 

(4.19 

-12 = y2i ■ ''S"- 

(4.20 

*21 " >'l2 * 

(4.2i; 

*22 ' >'22 ■ ■" ^S^y' 

(4.22 

*11 ' * Ip 

(4.23: 

*12 " *21 ' ''Vy 

(4.24 

^22 “ ^'^y^ 

(4.25 
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4.2. Temporal Discretization 

Following the procedure of Beam and Warming [4], the momentum 
equation (4.1) is now discretized in time. Using the superscript ”n‘' 
to denote a particular time step, the time increment of q is 


, n n+1 n 
Aq « q - q 


(4.27) 


The 3-point backward-difference expression for 3q/9t is 


. n+1 - n+1 , n . n-1 

(|f) '' 

Incorporation of Eqs. (4.27) and (4.28) into (4.1) leads to 




+ I At{^[Wq^ + Xq^]"'^^ + ^[Yqg + (4.29) 


+ 0(At3) 


Since A(q) and B(q) are independent of q_ and q , it follows that 

fe n 


} ■ ■ ■ 

aA""-^ 

3a’ 

r 

H 

" 9? 


aB”"^^ 

afi' 

■i 

3n 

9n 


where 


+ Q^Aq" + 0(At2) 

n 


* . / ^ n+lv 

A = A(q ; p ) 

n n+lv 
B = B(q ; p ) 


(4.3 

(4.3 

(4.3 

(4.3 
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Pij(q) 


P2i<q) P22(‘l> 


qil(q) qi2(q) 


q2i(q) q22(‘i> 


“ l(2?^u + y) 


6 u 

Pio(q) “ 


5 V 

poi (q) " “f" 


P22(q) = + 2y) 


qil(q) = 3-(2n^u + nyV) 


n u 

qioCq) =‘ 


n V 

■•21 ‘O) -i- 


q22(q) = + 2nyV) 


Now, assuming that p,p is a slowly varying function of time, it is 


admissible to replace Eq. (4.29) by 




a. ^ A.- 9 fv^A a , _n. n ^n. nT 
+ j At -^[Y Aq^^ Aq^ - Q Aq J 


. 2 . .f 3 r.,a n . „n n **t . 3 r^n n , „n n _*t-. 

+ 3 At{^CW q^ + X q^ - q^ + Z q^ - B ]} 


+ 0(At3) 





Following Beam and Warming, the cross derivatives are time-lagged 
in Eq. (4. 44) using the substitutions 


^<Xiq") . + 0(4t2) 


(4.45) 


IfWq") - ^OtAq"-*) + 0(4t^) 


(4.46) 


The other derivatives of Aq”, Aq”, and Aq** are moved to the left side 

4 n 

of (4.44), resulting in 


I At -|r(WAq" - P%") - f At -^(Z4q" - Q^q") 


iM 

3 J 


2 8 > n * 


(4.47) 


+ I At + Zq” - B*) + 0<^t3) 


where 


* n , . n-1 
q = q + Aq 


(4.48) 


4.3. Approximate Spatial Factorization 

After multiplying through by J, the left side of (4.47) is factored 
into a product of 5- and f|“ operations , leaving the following sequence 
of equations to be solved : 


^ + D^((P” - W ■^)Aq) = I Aq"'^ + D^(Wq" + Xq* 
+ D^(Yq^ + Zq“ - B ) 


(4.49) 


where 


Aq” + D,((q” - Z ^)Aq") = Aq + 0(At3) 

Ti on 


2 TA4- X 

3 9g 

2 ta.- 9 


(4.50) 


(4.51) 


(4.52) 




^ Jh. si:. ^ 


4.4. Spatial Discretization 

Spatial discretization o£ Eqs. (4.49) and (4.50) is accomplished 
by means of the central-difference expressions given in Appendix A. 
After some regrouping of terms, the spatial difference equations rep- 


resenting (4.49) and (4.50) take the form: 








“ 3 ^^ij^ "^^i+l,j ^ Vl,3^“^®i,j+1 " ®l,j-l) 


j-1 ^1, j+P^^i, j-1 ‘ ^^^i, j-1 ^i, j+l^‘‘ij ^^^i, j+1 ■*■ j-P**!, j+l 

^ k k ^11 k k 

+ - 9^+1 j_i) - 

* ~n * * 

■'■ ^i,j+l^%+l,j+l ■ ‘*1-1, j+1^ " \,j-l^‘*i+l,j-l “ 


and 


(4.54) 
‘ij 


■^^^i,j+l " ^^i,j+l ■ ^i,j-P^‘*l,j+l “ 


where 


p*" = 

At 

j , 


3 

ij 

In 

At 


Q * 

3 

"ij 


,n 


(4.55) 

(4.56) 


t The subscripts "i” and ”j" are indices for finite-difference mesh 
points, with 5=i, n” j , and AC=An=l. 
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w“ “ X ■'ij"" 

i“ - f 

in . ^ 


«.57) 

(4.58) 

(4.59) 

(4.60) 


I = 


1 0 
0 1 


(4.61) 


Note that in Eqs. (4.55) - (4.60) the Indices of are held fixed; 

e*8» > 


p*' a At J p*' 

^i+l,j 3 ^ij^i+l,j 


(4.62) 


4.5. ADI Sequence 

Equations (4.53) and (4.54) in sequence resemble one iteration 
of an alternating-direction implicit (ADI) scheme. First Eq. (4.53) is 
solved on each line of constant h (constant j); then Eq. (4.54) is 
solved on each line of constant g (constant 1). Actually, Eq. (4.53) 
represents a set of linear equations for the values of 4q^^, where j is 
specified and i varies from 1 to I^^^. The same applies to Eq. (4.54) 
for 4q^j , but with i specified and j varying. In both cases the matrix 
operators occur in block-tridiagonal form, and the solutions can be ob- 
tained directly (i.e., without iteration) using the well-known Thomas 
algorithm [8,12]. For a known pressure distribution, only one ADI iter- 
ation is necessary to obtain the velocity distribution, correct to 


second order in time. 


CHAPTER V 




I 


FORMULATION AND SOLUTION OP THE PRESSURE EQUATION 
5.1» Conservation of Mass 

There often arise situations in which volume changes are neglig- 
ible compared with other phenomena. In particular, liquids and gases 
flowing at low Mach numbers exhibit essentially the same behavior as 
Incompressible flow. 

From a mathematical standpoint, incompressibility eliminates the 
equation of state, which relates pressure, volume, and temperature. As 
a result, the thermal and kinetic energies are decoupled, and the flow 
can be analyzed on a purely mechanical basis, without reference to 
thermodynamic effects. 

Unfortunately, the incompressible equations of motion are somewhat 
Inconvenient for computational solution because nowhere does there ap- 
pear a time derivative of the pressure. Since there is no direct way 
of advancing the pressure in time, an indirect method must be formulated 
such that conservation of mass is always satisfied. In this regard, the 
continuity equation is a constraint that determines the instantaneous 
distribution of pressure throughout the flow field* 

5.2. Formulation of the Pressure Equation 

■ 

Taking the divergence of the momentum equation in certeslan coor- 
dinates leads to 


t The x-derivative of (2.17) plus the y-derivative of (2.18). 
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d + p + p + + 2u V 

t *^xx ■ yy X y y x 


*• rT“[M (u + u ) + M (v + V ) + M u 
Re‘’'^x' XX yy y xx yy xx x 


(5.1) 


+ n (u 4* V ) + M V ] 

'^xy y x' yy y 


where 

6 « u + V 5 0 (5.2) 

X y 

In obtaining Eq. (5.1), it is assumed that 5“0, but that The 6^. 

is a correction term (Ideally zero) recommended by Hirt and Harlow [13]. 
Note that Eq. (5,1) has been developed in nonconservative form, and d 
has been extracted and set equal to zero. Following Roache [ll], this 
can be taken one step further: 

u^ + v^ *» d^ - 2u V (5.3) 

X y X y 


Combining Eqs. (5.2) and (5.3) with (5.1), the pressure equation in 
cartesian coordinates then becomes 


p + p = - d. 4- 2(u V - u V ) 
^xx *^yy t X y y x 


+ (u 4- u ) + p (v 4- V ) 
Re‘'^x' XX yy y xx yy 


(5.4) 


4-u u 4-p (u 4-v)4-p v] 
*^xx X xy y x yy y 


Equation (5.4) is an elliptic partial differential equation (PDE) 
that establishes the instantaneous relation between pressure and veloc- 
ity everywhere in the fluid. Representing a linear combination of Eqs. 
(2.17) - (2. 19) , it is a restatement of conservation of mass, contin- 
gent upon conservation of momentum. Ideally it should be solved simul- 
taneously with Eqs. (2.17) and (2.18), but the approach taken herein is 
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to soXve the momentum and pressure equations separately, The vector 
momentum equation is of the parabolic type, and the approximate- 
factorizstion algorithm described in Chapter IV is currently the most 
efficient technique available for parabolic-hyperbolic PDE's. On the 
other hand, the pressure equation is elliptic and is unsuited to the 
method of Chapter IV because it lacks a time derivative, 

5.3, The Pressure Equation in Curvilinear Coordinates 

Since Eq. (5.4) is an auxiliary constraint and not a conservation 
equation per se, some leeway is admissible concerning its transformation 
to curvilinear coordinates. Accordingly, the left-hand side is ex« 
pressed in geometrically conservative form (section 3.2.), but the 
right-hand side is expressed in geometrically nonc onservative form, 
using the chain rule. The resulting equation is: 


where 


= ^ 5^ + 2(u V - u V ) 
t X y y X 

2 

‘ TT”[u (u + U ) + p (V + V _) 

Re'" X XX yy y xx yy 

''x’ ^ 

(5.5) 


(5.6) 


(5.7) 

z- j(t)2 + np 

(5.8) 


and the derivatives on the right-hand side of (5,5) are evaluated as 
follows 
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f * e f - + n f 
X ^x*4 ^x n 


(5.9) 


f g f - + n f 
y y ^ y n 


(5.10) 


* \«x>, 


(5.11) 


(5.13) 


'xy " iC«y”x>5 "y^x^ «x«y>5 + 

'yy " «y'V? + ^V„ <=•»= 

Note that Eq. (5. 12) is written such that f » f . 

xy yx 

5*4. Temporal Discretization 

Using a two-point backward-difference approximation for the time 


derivative, 


,n+l - 6" . 

^t ” If ■ - -t om 


(5.14) 


the desired velocity field (u,v) at time level "n+1" must satisfy con- 
tinuity, i.e. , 


5”’^^ = 0 


(5.15) 


Due to numerical error, however, 5“ is never precisely zero. Thus the 
two-point backward-time approxiroation for Eq. (5.5) is 


t Strictly speaking, a three-point (second-order) expression should 
be used to restrict the truncation error to O(At^), as in Chapter IV. 
This, in fact, was the original approach. After some computer experi- 
ments, however, it was found that the first-order pressure equation is 
much more stable, making truncation^error considerations academic. 




+ ^[Spf 1 + ip;;"^]) 


+ St''x<V + “yy) •*• '‘y(\x + ''yy> 


& \x u ^ u ( u + V ) + U V 1 
XX X ^xy' y x' ^yy y** 


n+i 


+ 0(AU) 


n+1 

In order to solve Eq» (5.16) for p i the velocity distribution 

ii+1 

(upv) must be known. Noting that 


n+1 n , A/A..S 
u " u + 0(At) 


(5.17) 




m v” i 0(At) 


(5,18) 


The velocity distribution (u,v)” can be substituted for (upv)'^^^ in 
Eq. (5.16), preserving the instantaneous truncation error 0(At). 

5.5. Spatial Discretization 

As in Chapter IV, the viscosity n is regarded as a slowly varying 

n+1 n 

function of time, so that p ** p . Using central-difference approx- 
imations for the spatial derivatives (Appendix A), the difference 
equation representing Eq, (5.16) is 


^ iiHhl ^ n4*l ^ tv^l ^ nHhl ^ n4’l 
A44P4 1 4 + B44P4T + • + 1>4 + E. .P"^*_^, 

ij l-l,j ij ij ij^^l+l.j 1^1,1 -1 ij‘i,j+l 

. * n+1 , „ n+1 * n+1 , (■ n+1 _ 

®ljPi+l,j+l %jPl-l,j+l ^ijPi-I,j-l “ ^ij 


(5.19) 


where 


"id 


^ (3W^_j, j -t 


( 5 . 20 ) 
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(5.21) 


ij 


J4 4(W. , 4 + W. , -f Z ) 

ij' x-l,j i+X,j l,j-i l,j+l 


C 44 • T J44<3W,^, , + W4_, 4 ) 


"Ij A "1-1, j 


1 


^44 “ '^4 4(3Z, , , + Z, 4 , , ) 


Ij 4 “i,j+l' 


1 


E 44 “ T J 44 OZ 4 4 ^, + Z 4 4 _,) 


ij 4 ij"'“i,d+i ’ “i,j-r 


1 


®4 4 ~ /. 34 4 (X 4 i, 4 + Xi 4 . 1 ) 


Ij 4 “ij'“i+l,j '‘i,j+r 


\j “ ■ 4 '^ij^Vl.j \,j+P 


“"4 ^i,j-l^ 


1-4 •“ T J 44 <X 4 , . + X 4 4 ,) 

ij 4 ij i-l,j i,j-l 


(5.22) 

(5.23) 

(5.24) 

(5.25) 

(5.26) 

(5.27) 

(5.28) 


and 


^n 


in Ij . 0 / rv n Hv 

F , . = - 7 -^ + 2 (u V - u V ) 
ij At ' X y y x ij 


+ ‘fT’[p'^(u'^ + u" ) + y’^(v" + ) 

P.e*’ X XX yy y' xx yy 


(5.29) 


. n n , n / n , nv . n nn 

+ u u + V (u + V ) + y V J 
^xx X "^xy y x' yy y'ij 


The X- and y-derivatives in Eq. (5.29) are evaluated using Eqs. (5.9) - 
(5.13), and the and n~ derivatives appearing therein are calculated 
numerically from Eq. (A. 1) . 


t Note that the quantity 6 ^^ represents here the divergence of velocity 


evaluated at |- 1 , n= j > and is not to be confused with the Kronecker 
delta used in Eq. (2.1). 
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5 , 6 , Iterative Solution 


Equation (5.19) reptosents ^ of N linear equations in N un- 


n+1. 


knoOTis ), where N is the total number of finite-difference node 


points in the flow field. For large N, direct solution is impractical 
and inaccurate, leaving indirect (iterative) methods as the more viable 
class of alternatives. Of the latter, successive over-relaxation by 
lines (line SOR) has been chosen as the method of preference, because 
it converges rapidly for elliptic equations with both Neumann and 
Dirlchlet boundary conditions (see Chapter VI). 


Dropping the superscript "n-fl" and replacing it with "(m)", to de- 
th 


note the in— iteration, the following equation is solved on each line 
of constant 5 (constant i): 


A ^ A ^ A ^ 

°ljPi,j-l ” ®ijPij " \jPi,j+l 


r X A (l«-l) ^ n (m-1) 


(5.30) 


. : (m-i) (m-i) 

°id‘^i+i,j+i ^ Pi-1, j+i 


+ K + r 


Equation (5.30) actually represents a set of linear equations for the 


value of p^j at each value of j . The matrix operator for the lef t-hand 


side is point tridiagonal, so that Eq . (5.30) can be solved directly 


using the Thomas algorithm. Once the values of p^^ have been found on 


(m) . 


a given C-line, the values of p^^ are calculated from 
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( 5 . 31 ) 

where the acceleration parameter m must lie between 0 and 2 for con- 
vergence. 

The foregoing procedure is repeated for each |-line in succes- 
sion, until the entire field has been swept. One sweep of the field 
constitutes a single iteration of line SOR. 




CHAPTER VI 


INITIAI, VALUES AND BOUNDARY CONDITIONS 

6.1. General Requirements 

In order to obtain a unique solution for a system of partial dif- 
ferential equations, it is necessaty to specify certain initial and 
boundary values for the unknown functions and their derivatives, de- 
pending on the type(s) of equations involved. In particular* for the 
incompressible Navier-Stokes equations, the Initial velocity and pres- 
sure distributions must be given, and the velocity (or a spatial de- 
rivative thereof) must be specified on the physical boundaries at all 
times . 

The no-slip condition for viscous flow sets the relative velocity 
at 2 ero on all fixed impermeable boundaries (in this case the airfoil 
surface). Freestream velocities may be specified at will, subject to 
conservation of mass and momentum. 

There is no self-evident boundary condition for the pressure, nor 
is one required in the strict sense. In fact, the pressure can be elim- 
inated entirely by using a stream-function/vorticity formulation in two 
dimensions [ll]. Nevertheless, when the primitive variables (u,v,p) 
are retained In the Navier-Stokes equations, and a Poisson equation is 
to be solved for the pressure, then a pressure boundary condition is re- 
quired. 

6.2. Freestream Boundary 

Ideally, the freestream (inflow/outflow) boundary lies an infinite 
distance from any obstacles in the flow. On a computational grid, 

"infinite" means far enough away that their effect on the flow is weak. 
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Here a freestream velocity ox a freestream pressure gradient can be 
specified (one Implies the other due to conservation of momentum) and 
the flow Is usually presumed to be Inviscld. When the upstream con- 
dition is that of uniform flow, it is customary to hold the velocity 
and pressure constant on the inflow boundary. On the outflow boundary, 
either the flow variables or their gradients may be specified. 

6*3. Gradual Start 

The inflow is accelerated from zero to its final velocity 
(|u|^ *= 1) by imposing a uniform body force on the entire flow 
field. During this phase, the Inflow velocity is given by 


u s /g dt 

•ycci ° 


( 6 . 1 ) 


and the outflow boundary conditioiT is 

(n*Vu)^ =0 (6.2) 

where n is a vector normal to the outflow boundary. The freestream 

fw 00 

pressure condition is 


Vp^= 0 (6.3) 

and g is a vector of uniform magnitude and direction, so that 

y-g = 0 (6.4) 

Thus the pressure equation (5. 5) is unaltered by the presence of g . 
t^hen the acceleration phase is complete, the body force is of course 

g =0 (6.5) 
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Since the freestream boundary is far from the airfoil, Eq. (6.3) can be 
replaced by the condition 

* 0 ( 6 . 6 ) 

The inflow velocity condition after the gradual start is 

u cos 01 (6.7) 

V ■= sin a (6.8) 

where a is the angle of attack. 

The body force is equivalent to a freestream pressure gradient 
acting on the fluid, and it can be treated as such, both analytically 
and computationally. Thus the flow could be accelerated by adding a 
static pressure gradient to the right-hand side of the momentum 
equation. In any case, the body force influences the field pressure 
through the pressure boundary condition on the airfoil surface, regard- 
less of the interpretation. 

6.4. Pressure Boundary Conditio n 

The no-slip condition requires that u = 0 on the body surface, so 
that the normal component of the momentum equation (including the body 
force) reduces to the nonconservative form 

n*Vp = n* (g + ~ V^u) (6.9) 

This equation is Simplified further by the boundary-layer approximation, 

n*Vp = n*g (6.10) 
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which eliminates the coupling between the pressure boundary condition 
and the velocity field. 

In curvilinear coordinates, the airfoil surface is a line of con- 
stant n» and 7n is a vector normal thereto. Thus Eq. (6.10) is equiva- 
lent to 

VirVp “ 8‘Vn (6.11) 

fV . 

Using the chain rule to expand 7ri and Vp, Eq. (6.11) becomes 

where gj^ and 82 are the x- and y-components of g. 

Equation (6.12) constitutes a Neumann boundary condition for the 
pressure on the body surface. In order for the llne-SOR calculation of 
the pressure to converge, Eq. (6.12) must be incorporated into the dif- 
fetence equation (5.19) evaluated either on the body or one line off 
the body. As Roache [ll] emphatically points out, it is insufficient 
simply to extrapolate the body values from the field without first mod- 
ifying the coefficients in Eq. (5.19), subject to Eq. ( 6 . 12). The most 
straightforward approach, eliminating the need for extrapolation alto- 
gether, is to evaluate and solve the pressure equation on the body 
simultaneously with that in the field. ^ 


+ All attempts to use Eq. (6.9) as the pressure boundary condition 
led to divergent pressure solutions, apparently because the pres- 
sure and velocity calculations are not done simultaneously. 
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Letting the n~line for j“s represent the airfoil surface, and 
using central-difference approximations (Appendix A) for p^ and p^ in 
Eq. (6.12), the following expression is obtained for p^ 


(? n + C n ) 

X X y y 


Pl,s-l “ Pi,s+1 ( 2 + J) ■‘’i+l.s " Pi-1, s^ 

' X y' 


(6.13) 


- ^^Vl V2^ 

X y 


Using first-order one-sided difference approximations (Appendix A) for 
the n-derivatives in Eq. (5.5), and combining the result with Eq. 
(6.13), the difference equation for the pressure on the body reduces to 


where 


t 0+1 , ^ n+1 . n+1 

^1,3 Pi-1, s l,s Pl,s i,s Pi+l,s 

, ^ n+1 , p „o+l . u - w 

i,s Pi,s+1 i,s Pi+l,s+l i,s Pi-1, s+1 i,s 


(6.14) 


K.S - + I *i.S> 

2(n^ + n^) 


(6.15) 


Bj “ - (W, , , + W. , + 2Z, 

i,s ' i+l,s 1-1, s l,s+l 


(6.16) 


4 - / , J - 4-(X, ., ^ + X, ■ ) 

i,s 4 i+l,s i-l,s 2 i+l,s i,s 




+ 

X 'y 


h + Z. .,) 

i,s i,s+l' 


(6.17) 


E, = 2Z. . , 

i,s i,s+l 


(6.18) 
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- K+I.S + *l,s+l> 


(6.19) 


and 


H 


1,8 


■ ■*■ \,s+P 


(n^ Si ^ * 

l.S / 2 4. .2) i*s i,s+l' 

\ ly/ 


( 6 . 20 ) 


( 6 . 21 ) 


Note that 6^ has been set equal to zeto on the body. 

6.5. Re-entrant Boundaries 

Re-entrant boundaries occur wherever a branch cut is made In the 
physical plane (Figure 1). Velocity and pressure and gradients there- 
of must be continuous across these boundaries, so the re-entrant bound- 
ary conditions are periodic in the transfomjed plane. To avoid extrap- 
olating these boundary values from the field, the physical variables on 
both sides of the cut are calculated simultaneously. That is, two 
5-llnes that meet on the cut represent a single ^-line as far as the 
solution algorithms are concerned, and the q-dlrection matrix inversion 
Is continuous across the wake. 

6.6. Trailing Edge 

The body-fitted coordinate system shown in Figure 1 is of the 

C-type or wake type, in which the iT'llaes begin and end on the outflow 

+ 

boundary. If the airfoil Surface coincides with the line p-s, then the 
two points and f;=r occur at the trailing edge, and the re-entrant 
segments are defined by F, ^ i and 5 r along n=s. 


t As Opposed to an 0-typc coordinate system, in which the 
p-lines are closed curves. 
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Ac Che Crailing edge, which is a sharp poinC in Che physical 
plane, Che surface-normal vecCor (Vn) is disconcinuous , making Che 
pressure boundary condition (6.11) also disconcinuous. This singular- 
icy is peculiar Co Che C-Cype coordinate sysCem rather Chan Che flow 
iCsclf, and it can be clrcumvenCcd on physical grounds. 

The no-slip condicion implies Chat cherc can be no unbalanced 
force on Che fluid aC Che Crailing edge. NeglecCing Che viscous sCress, 
as in secCion 6,4., Che projecCion of Che momenCum equaCion in an ar- 
biCrary direcCion t reduces Co 

T*7p « T*§ (6.22) 

If t is Cangent Co Che re-enCranC segmenCs aC and C»r on n*s, then 
the pressure boundary condition aC the Crailing edge becomes 

P^ * j”Unygj - njjg2> (6.23) 

as and ^‘♦■r, . 

“■ *r 

CompuCaCionally, the following substitutions are made in the 
line-SOR equation (5.30) at the trailing edge. 


'i+l,s " Pjl-l,s “ \®2^ 

(6.24) 

r-l,s * **rfl,s ~ T^*^y®l ~ \®2^ 

(6.25) 
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CHAPTER VII 


computational adjustments 

% 

7.1. Turbulence Model 

Baldwin and Lomax [14] have formulated an algebraic turbulence 
model for separated flow, which is used herein to Include the effects 
of small-scale eddies in the Reynolds-averaged Navier-Stokes equations. 
This model generates an eddy viscosity ^ 0 such that the nondimen- 
sional physical viscosity becomes 

M = 1 + Prp (7.1) 

and *>■ 0 on the body and in the freestream. Transition is set at the 
point of minimum pressure on either side of the. stagnation point, with 
the turbulence propagating downstream into the wake. 

7.2, Artificial Viscosity 

Space-centered differencing leads to algorithms that are algebraic 
ally convenient, second-order accurate, and unfortunately susceptible 
to Instability, especially at high Reynolds number. Some form of arti- 
ficial dissipation is usually needed to diffuse spurious oscillations, 
even with fully implicit methods, This task is accomplished herein by 
adding an extra viscous term to the (nondlmenslonal) right-hand side of 
Eq. (2.12); specifically^, 

WS(2.12) = RHS(2,12) + + jjl) (7.2) 

j i 1 

The additional term is the same as the viscous term in (2.12), but with 


t Equation (7.2) uses the Einstein summation convention and the 
cartesian coordinates of Chapter II. 
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ji“l and Re replaced b’j? the artificial Reynolds number [15], 

defined by , , 

e|v-u| 

Ra “ J 


(7.3) 


where e ^ 0 and J is the Jacobian of the coordinate transformation, 
given by Eqs, (3.16) and (3.18). 

Ideally, the artificial viscous term should vanish wherever the 
calculated flow satisfies continuity. The main advantage of (7.3), 
however, is that it behaves like a switched filter, turning on wherever 
instabilities are most likely to arise. For situations in which *7*u 
varies sharply across the field, an extended definition of Ra is 



where the subscript "ave" indicates a three-point average in the di- 
rection of sharp variation. In either case, the additional term in 
(7.2) introduces both implicit and explicit artificial viscosity into 
the approximate-factorization scheme discussed in Chapter IV. 

7.3. Pressure Smoother 

The primary source of instability in the calculated flow is the 
pressure equation. Regardless of the method of solution or the degree 
of convergence, the right-hand side of the difference equation (5.19) 
contains only information from the previous time step. As a result, 
numerical oscillations In the velocity field are fed directly into the 
pressure equation, worsening the situation with each time step. 

Viscosity, either real or artificial, dissipates oscillations in 
both pressure and velocity. The need for artificial viscosity can, how- 
ever, be reduced somewhat by smoothing the right-hand side of the 
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pressure equation. Moreover, a three^polnt average In the direction o£ 
sharp variation improves results considerably. Specifically, Eq. 

(5.19) is replaced by 

UlS(5.19)-{(Pj^,_j+2»y+F^.j_j) (7.5) 

which is more or less equivalent to the addition of a term proportional 
to '?*‘*p in the pressure equation. 
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CHAPTER VIII 
COMPUTATIONAL RESULTS 

8.1, Coordinate Sy at em for NACA 66 ^ 018 Airfoil 

Thompaon [15] Has generated a body-fitted curvilinear coordinate 
system lor the NACA IlgOlS airfoil, based on an assumed Reynolds num- 
ber of 100,000, This coordinate system, shown in Figure 2, is suitoble 
for calculations at Reynolds numbers within an order of magnitude or so 
of 100,000, If Re » 100,000, the flow resolution in the viscous region 
will be poor, because the mesh spacing is not fine enough next to the 
body. On the other hand, if Re « 100,000, round-off problems may 
arise because the mesh spacing is too fine. The transformed plane is 
a 113 X 51 rectangular grid with unit mesh sp'ieing, 

8.2, Computational Parameters 

In all the calculations reported below, the time step is fixed at 
At«0,01. The gradual start consists of 100 time steps of uniform accel- 
eration, with the X and y body-force components given respectively by 

gj^ " cos a 

$2 “ sin a 

After the first 100 time steps, at t«1.00, the calculation is stopped 
and then restarted with the initial condition Au =■ Av ® 0, Since a 
three-point backward time difference is used in the AF momentum solver 
(section 4,2,), the testart is necessary to avoid overshooting the val- 
ues of Au and Av at -‘»l.0l. Had a nonlinear (e,g,, cosine) aoceleration 
been used, the stop and restart would be unnecessary, because the body 

force would be a continuous function of time. As it is, the body force 
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for the linear start is a step function, discontinuous at t"0*0 and 
t«*1.0. 

The number of iterations in the line-SOR pressure solver is 
fixed at thirty per time step, with the acceleration parameter set at 
w**1.85. Thirty iterations are probably more than necessary, although 
stability does improve slightly with the iteration count. Reducing the 
number to twenty degrades stability and convergence only slightly; 
further reduction to ten or fewer seems to promote instability, due to 
poor convergence. The oj-value of 1,85 is not necessarily optimal, but 
rather the result of some experimentation. Values from 1,8 to 0,0 
lead to increasingly slower convergence, while the optimum seems to lie 
somewhere between 1.8 and 1.9. 

The artificial viscosity coefficient can be sec at any value e> 0 , 
depending on how stable the undamped calculated flow is. Since this 
parameter alters the local Reynolds number (section 7.2.), it is de- 
sirable to keep the value of e as low as possible. The pressure 
smoother (section 7.3.) is used in all calculations, since it does not 
add any fictii/ous terms directly to the momentum equation, nor does it 
affect the spatial truncation error. 

8.3. Results for Re ~ 1000 and a = 6° 

Even though the coordinate system (Figure 2) was designed for 
Re =100,000, it is worthwhile to make a calculation at low Reynolds 
number, without including the turbulence model and the artificial vis- 
eosity. The object here is to examine the qualitative behavior of the 
predicted flow field in the absence of artificial stabilizing mechanisms. 

The results presented in Figures 3 - 12 were obtained for the con- 
ditions 
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Flgurea 3 6 show sequentially the ovolvielon oC the Rreasuee til8ti?l))n- 

tlon along the alffollt^ I)ownsteoaiu ol tlto loading edge, tho pr^esiu*© 
on the lower (windward) surface gradually falls, while, ulto presaure on 
the upper surface gradually rises. At t“3.0 the two distrihutiona have 
begun to ovorlap, with the point of intersection moving upstream to 


x'’0«5 by t*4.0. 

Figures 7 “ 12 Illustrate the developmont of tlie velocity field 
between t“2.0 and t*4.0. Around the leading edge (Figures 7 and 8), 
there is hardly any change in the velocity profiles. Moreover, exami- 
nation of tho entire flow field (Figures 9 and .10) reveals little 


change except near the trailing edge CFigures 11 and 12). Hero, a vor- 
tex Ju.st beginning to form at t*2.0 has grovvn considerably by t*»4.0. 

The calculation was tciMnated at t«4.3, duo to Instability asso- 
ciated w4tb the trailing-edgo vortex. A restart could bavo been initi- 
ated at t*4,0, with artificial viscosity acklod for stabilization, but 
this seemed unv«irran ted. The informatiou generated for t<4*Q gives an 
adequate picture of tlio flow at low Reynolcls number; and the instabil- 
ity Is probably due to the meBlr s]>aeing, which is imieh too fine. 


1‘ The loading edge coincides with x**0, and the trailing edge with x®l. 
The pressure cooffieient C is defined by 


p 


C 


P 


1 _ it2 
'o p i\s 




8*4. Results for Re ° 40,000 and a ” 0 


Mueller [16] has conducted wind-tunnel tests of the NACA 
airfoil at Reynolds numbers from 40,000 to 400,000. In order to eval- 
uate the applicability of the combined AF/SOR scheme for real flow 
problems, calculations have been made for comparison with Mueller's 
data. The computational results presented in Figures 13 - 34 were ob- 
tained for the conditions 

Re '= 40,000 
e « 1 , t S 1 • 0 

E =* 10, t > 1.0 
a = 0° and 6° 

Trial calculations with e=0 became unstable around the leading 
edge at t=0.8. With E=l, stability problems developed in the wake, 
near the trailing edge, at t~1.8. Thus, to maintain stability for 
Re = 40,000, it is necessary to increase the artificial viscosity coef- 
ficient from 1 Co 10 when t>l. Furthermore, to keep the flow 
well-behaved in the wake, it is advisable to activate the turbulence 
model prior to t=3.0.*t‘ 

The computation for a=0 CFigures 13 - 23) was executed from t=0.0 
to t=2.0 with = 0. The turbulence model was activated at t=2.01, 
and execution continued until t=5.0 (a total of 500 time steps). The 
combined eddy/artiflcial viscosity kept the flow orderly, and no 


t Since the turbulence model increases the local viscosity, it has a 

stabilizing effect on the flow, especially in the wake. As a result, 

the turbulence helps to suppress the development of unstable trailing- 

edge vortices of the kind discussed in section 8. 2.- (Figure 12). 
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instabilities were apparent when the calculation was finally stopped. 

Figures 13 - 17 indicate the change in the pressure distribution 
between u^l.O and t=5.0. The calculated variation of is syimnetric 
along the airfoil (as it should be), and the pressure drop at t®3.0, 
between x«0.8 and X“1.0, is precipitated by vortices at the trailing 
edge (Figure 22). These vortices have begun to damp out by t«5.0. 

The agreement between the computed and the wind-tunnel data is 
good everywhere except near the leading and trailing edges. Further 
calculation would probably eliminate the discrepancy at the trailing 
edge, since the computed pressure seems to be falling there with time. 
The leading-edge problem may be the indirect result of continuity vio- 
lations and artificial viscosity, via Eq, (7.4) , and might or might 
not be corrected by additional computation. 

Figures 18 - 23 show the change in the velocity field between 
t=2.0 and t=5.0. For the most part, the viscous layer around the body 
is unr< 2 alistically thick, with separation occurring at center chord 
(Figure 19). This is unavoidable with the AF/SOR method in its pres- 
ent form, since it is caused by the artificial viscosity, which is nec- 
essary for stability at high Reynolds number. 

At t«2.0 the leading-edge velocity profiles (Figure 20) are 
well-behaved, tbough getting thicker in the x-direction. By t=5.0., 
the profile shapes (Figure 21) have begun to alternate somewhat, due to 
the variation of the artificial Viscosity along the body. This 

t Experimentation with the AF/SOR method has indicated that the greater 

the penetration (continuity violation) around the leading edge, the 

gentler the pressure gradient there. 
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alternation could be the harbinger of trouble for t>5,0, but the 
velocity and pressure distributions are still synunetrlc and quite 
stable at t **5.0, 

8.5. Results for Re 40,000 and a ■* 6” 

The computation for a “6° (Figures 24 - 34) was executed from 
t^O.O until t«5,0, with the turbulence model activated in the second 
time step (t®0.02). No instabilities were encountered, nor were any 
apparent when the calculation was terminated. 

The evolution of the pressure distribution along the body 
(Figures 24 - 28) is qualitatively the same as for Re = 1000. In this 
case, however, the slopes are steeper, and the overlap for the upper 
and lower surfaces takes longer to develop. Comparison with the 
v^ind^tunnel data in Figure 28 shows the computed curves to be mysteri- 
ously displaced from the experimental curves, but with the same general 
Shape, Including the double overlap. Moreover, in Figures 28 and 17 
alike, the calculated pressure distributions appear to be right-shifted 
from the experimental results. In the 6-dcgree case, the shift may be 
due to poor resolution of the calculated minimum pressure on the upper 
surface near the leading edge. 

Examination of the velocity field (Figures 29 and 30) reveals 
substantial change between t=2.0 and t = 5.0. The upper- and 
lower-surface separation points, located respectively at x« 0,6 and 
x“ 0*8 when t= 2.0, have moved upstream to x - 0,3 and x« 0.7 by t = 5.0. 
The leading-edge velocity profiles exhibit a gradual thickening 
(Figuras 31 and 32), and vortices have developed at the trailing edge 
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by t»5.0 (Figures 33 and 34) • As in the calculation for a“0, the 
viscous loycr is much too thick around the body. 

8.6. Other Calculations 

Prior to making the calculations for the airfoil, a great 

deal of personal effort and computer time was Invested in calculations 
for the NACA 64A010 airfoil at Re “ 2,000,000, using a coordinate 
system generated by Cooper [S], In fact, most of the developmental 

4 * 

work for the AF/SOR scheme was done with the 64A010 airfoil. ' To main- 
tain stability with this coordinate system and Reynolds number, how- 
ever, the required values of & were 10 during the gradual start and 100 
thereafter. In other words, the amount of artificial viscosity needed 
for Re » 2,000,000 was an order of magnitude greater than that for 
40,000, making the viscous region about three times thicker at the 
higher Reynolds number. 

The stability problem is probably related to the coordinate geom- 
etry, specifically to the mesh spacing. Equation (3.11) determines the 
n-llne spacing next to the body; so if Re increases by two orders of 
magnitude, the spacing decreases by one order of magnitude. It is com- 
mon knowledge that, for explicit methods, the maximum stable time step 
decreases with the minimum finitc-diff erence cell dimension. The same 


t The AF/SOR method, as it now stands, represents the last of three at- 
tempts to construct an approximate-factorization algorithm for incom- 
pressible flow. The first two versions, while invaluable from an ex- 
periential standpoint, produced little in the way of communicable re- 


sults. 


is apparently true for Implicit/expllcit methods such as AF/SOR, but 
the relationship is more difficult to quantify. In either case, it 
seems that more viscosity may be needed for stability when the mesh 
spacing is reduced with the time step and Reynolds number fixed. Ad 
ditional calculations are needed to verify this for AF/SOR. 
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CHAPTER IX 
CONCLUSION 


9.1. Discussion 

The combined AF/SOR algorltlun Is Implicit In its solution of the 
momentum equation, but explicit In the sense that the pressure/ 
continuity solution is lagged in time. As a result, the method is far 
less stable than the compressible-flow AF scheme of Beam and Warming. 
For a time step At = O.Ol, the maximum practical Reynolds number is 
probably around Re =10,000, If the artificial viscosity is to be kept 
tolerably low (s ± 1). Furthermore, since stability may depend on 
mesh spacing as well as viscosity, the body-fitted coordinate system 
should be generated for the Reynolds number in question or for a lower 
value, but not a higher. Increasing the Re-value in the coordinate gen 
eration (section 3.1.) reduces the mesh spacing. This, in turn, de- 
mands more artificial viscosity for stability, which contaminates the 
numerical results. Since truncation error is preferable to instabil- 
ity, it seems advisable to keep the mesh a little on the coarse side. 

By way of comparison with conventional implicit methods, point-SOR 
calculations can be stabilized at Re - 10® with e = 1 and At = 0.01. 

Thus, at high Reynolds number, point SOR (or some other fully Implicit 
method) is required for incompressible flow calculations. On the other 
hand, at low Reynolds number (Re < 10**) AF/SOR seems preferable, being 
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passably stable ond fastec than point SOR by at least a factor of tvo.^ 

Compared with compressible-flow AF calculations, AF/SOR fairs 
poorly. Cooper [8] reports a Machine-time expenditure of 4.5 seconds 
per time step on the CDC Cyber 203, using the method of Beam and 
Warming at a freestream Mach number of 0.8 and a Reynolds number of 
two million. The AF/SOR calculations discussed in Chapter VIII aver- 
aged 5.5 seconds per time step on the same computer. While program 
optimization could markedly reduce both of these cost figures, the 
point is that the compressible-flow algorithm seems to require no more 
computer time than the Incompressible solution. 

In summary, the combination of approximate factorization and line 
SOR appears to be suitable only for incompressible flow calculations at 
low Reynolds number (Re < 10^*), when used alone. On the positive side, 
however, AF/SOR might serve as an initial-guess generator for iterative 
methods such as point SOR. It is possible that the latter approach, 
employed in each time step, would produce stable high-Reynolds-number 
solutions in relatively few iterations compared to point SOR alone. 


t Thirty iterations of line SOR is roughly equivalent to fifty itera- 
tions of point SOR, for a single equation. Ten iterations of the 
line-SOR pressure solver takes roughly the same number of operations as 
the AF momentum solver. The net result is that AF/SOR, as used herein, 
requires between two and three times less work than a 50-iteration 
point-SOR solution for u, v, and p. 


9.2. Epilogue 


The results presented in Chapter VIII are in some ways disap- 
pointing, but not really surprising. While fully explicit and fully 
implicit methods usually have definable (or at least estimable) limits 
of stability, mixed methods can only be evaluated by computational ex** 
perlmentation. At the outset, one hopes for the best, but in the end 
must accept those limitations dictated by the ^Iachlne• 

Incompressible flow at high Reynolds number remains a knotty prob- 
lem, requiring expensive, and sometimes questionable, methods of solu- 
tion. At low Reynolds number, AP/SOR offers a means of calculation 
that is relatively efficient. Unfortunately, most engineering problems 
involve high Reynolds number, whether they be aerodynamic or hydrody- 
namic. Only rarely are calculations needed for molasses or for concen- 
trated shampoo. 

What is needed is a method that calculates pressure and velocity 

simultaneously, but with the same efficiency as approximate factoriza- 

i* 

tion. A vectarized ADI or SOR scheme could meet the simultaneity cri- 
terion but would be inefficient, requiring far more than one iteration 
per time step for convergence. Perhaps a clever formulation or manipu- 
lation of the pressure/ continuity equation would allow factorization in 
the same fashion as compressible flow. Except for the method of Steger 

t That is, an algorithm that inverts a 3 x 3 matrix at each point in 
the field to obtain the vector (u,v,p) . This takes about four times 
as many operations as the more commonly used approach, which is to cal- 
culate u, v, and p in sequence at each point. The latter method, in 

effect. Inverts three 1x1 matrices. 
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and Kutlar [9]» which converges only for the steady state, efforts in 
that direction have failed to produce a stable algorithm. Neverthe- 
less, the problem is so intriguing, so deceptively simple at a glance 
that it can become an obsession. Every day one thinks anew, there 
must be a way. 
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Figure 5. Pressure Distribution - Re « 1000, ct ' 


* 


57 



Figure 6. Pressure Distribution - Re - 1000, a = 6 
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Figure 7. Velocity Field - Re ■ 1000, a 


59 








nv' 


TO 


Figure 9. Leading --Edge Velocity Profiles - Re = 1000, a 





Figux«2 11 f Tralllng-Edge Veloctty Profiles - Re = 1000, a 





Figure 14. Pressure Distribution - Re « 40,000, a «= 0 
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Figure 15. Pressure Distribution 
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Figure 17. Pressure Distribution - Re = 40,000, a 0 
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Figure 21. Leading-Edge Velocity Profiles - Re = 40,000, 0=0°, 
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Figure 22. Trailing-Edge Velocity Profiles - Re = 40,000, o 
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Figure 24. Pressure Distribution ^ Re = 40,000, a 






Figure 25* Pressure Distribution - Re = 40,000, a 





Figure 27. Pressure Distribution 








Figure 29 . Velocity Field 













APPENDIX A 

FINITE-DiFFERENCE APPROXIMATIONS 


The difference expressions listed below were used In discretizing 
the spatial derivatives In the main text (except where otherwise 
noted). The grid spacing In all cases is AC “ An ® 1. Whenever the 
forms for the C- and n-derivatives are Identical, only one form is 
given. The Indices "1" and "j" denote C- and n”values, respectively. 
Second-Order Central Differences; 


“ 2 ^^ 1+1 " ^ 1 - 1 ^ 


(A.l) 




f - 2f + f 
1+1 1 1-1 


(A. 2) 




(A. 3) 

4^^i+i,j+i *■ ^1-1, j+1 " ^l+l,j-l ^1-1, j-l^ 

~ ^i-l,j+l^ "^i,j-l^^i+l,j-l “ ^1-1, 

(A.5) 

" i^+l,j^^i+l,j+l " ^i+i,j-P ■ Vl,j^^i-i,j+l " ^1-1, j-P^ 

(A. 6) 

First-Order One-Sided Differences; 
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(A. 7) 
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APPENDIX B 

A diffekencing problem 


Cooper [S] points out that the same type of differencing should 
be used in the flow calculation as in the numerical calculation of thn 
metric coefficients. For example, if central differences are used to 
calculate | and E » then the same should be used for Uj. and v„ , and 
so on. 

An additional problem has been encountered with respect to the 
viscous terms* Consider the two possible central-difference expres- 
sions for (Zu ) 


(Zu^) = 






(B.l) 




(Zu ) » (Zu. ) - (Zu ) . . . 

n fi n n 


^ 2^h,j+i ^i,j^“i,j+i 


1 , 
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(B.2) 


For some reason, Eq. (B*l) works and (B.2) does not. 'The same 

problem applies to (Wu») and, o£ course, to similar terms Involving 

b $ 

the other velocity component v. 

Wlien Eq, (B.l) Is used in the AF momentum solver, without over 
calling the pressure solver, the results ore stable and well*-bohaved 
(though erroneous next to the airfoil). On the other hand, if Eq. 

(B,2) is used, especially during the gradual start, oscillations arise 
in Che velocities. This hoppens with or without the pressure solver , 
and the situation .gets worse If the Reynolds number is reduced (just 
the opposite of what would be expected). 

The reason for the instability with Eq. (B.2) is not clear. The 
metric quantities n and n use values of x and y at j+1 and j-1 but 
not at j. It is possible that the coefficients in the difference ex- 
pression for (Zu ) are restricted in the same way. In any case, the 
n ,1 

remedy is to use Eq. (B.l) , not Eq. (B.2). 
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